arXi V: 1504.02891v 1 [math.NA] 11 Apr 2015 


A REGULARIZED NEWTON METHOD FOR COMPUTING 
GROUND STATES OF BOSE-EINSTEIN CONDENSATES* 

XINMING WUt, ZAIWEN WEN*, AND WEIZHU BAO§ 


Abstract. In this paper, we propose a regularized Newton method for computing ground states 
of Bose-Einstein condensates (BECs), which can be formulated as an energy minimization problem 
with a spherical constraint. The energy functional and constraint are discretized by either the fi¬ 
nite difference, or sine or Fourier pseudospectral discretization schemes and thus the original infinite 
dimensional nonconvex minimization problem is approximated by a finite dimensional constrained 
nonconvex minimization problem. Then an initial solution is first constructed by using a feasible gra¬ 
dient type method, which is an explicit scheme and maintains the spherical constraint automatically. 
To accelerate the convergence of the gradient type method, we approximate the energy functional 
by its second-order Taylor expansion with a regularized term at each Newton iteration and adopt a 
cascadic multigrid technique for selecting initial data. It leads to a standard trust-region subprob¬ 
lem and we solve it again by the feasible gradient type method. The convergence of the regularized 
Newton method is established by adjusting the regularization parameter as the standard trust-region 
strategy. Extensive numerical experiments on challenging examples, including a BEC in three di¬ 
mensions with an optical lattice potential and rotating BECs in two dimensions with rapid rotation 
and strongly repulsive interaction, show that our method is efficient, accurate and robust. 


Key words. Bose-Einstein condensation, Gross-Pitaevskii equation, ground state, energy func¬ 
tional, spherical constraint, gradient type method, regularized Newton method. 

1. Introduction. Since the first experimental realization in dilute bosonic atomic 
gases [3[22llSl], Bose-Einstein condensation (BEC) has attracted great interest in 
the atomic, molecule and optical (AMO) physics community and condense matter 
community [31133 iu sg. The properties of the condensate at zero or very low 
temperature are well described by the nonlinear Schrodinger equation (NLSE) for the 
macroscopic wave function 7 /^ = '0(x, t), which is also known as the Gross-Pitaevskii 
equation (GPE) in three dimensions (3D) [6l [29l [36l |43l HU [45] as 



( 1 . 1 ) 


where t is time, x = (x,^, z)^ G is the spatial coordinate vector, m is the atomic 
mass, h is the Planck constant, N is the number of atoms in the condensate, Q is 
an angular velocity, E(x) is an external trapping potential. The term Uq = 
describes the interaction between atoms in the condensate with the s-wave scattering 
length as (positive for repulsive interaction and negative for attractive interaction) 
and 


Lz = xpy - yPx = -ih{xdy - ydx) 
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is the 2 ;-component of the angular momentum L = x x P with the momentum oper¬ 
ator P = —ihV = {PxiPy^Pz)~^‘ If is also necessary to normalize the wave function 
properly, i.e., 

(1.2) ||^(.,t)||2 ;= / \^{x,t)\‘^dx = 1. 

By using a proper nondimensionalization and dimension reduction in some lim¬ 
iting trapping frequency regimes [191131], we can obtain the dimensionless GPE in 
d-dimensions (d = 1,2,3 when O = 0 for a non-rotating BEG and d = 2,3 when 
7 ^ 0 for a rotating BEG) [ini [43l[45] : 

+ E(x)-h/3|'0(x, I'0(x, t), X G t > 0, 


(1.3) 


dt 


with the normalization condition 

(1.4) ||V;(.,t)|p:= / |V^(x,t)pdx = l, 

jRd 


where /3 G M is the dimensionless interaction coefficient, Lz = —i{xdy — ydx) and 
V (x) is a dimensionless real-valued external trapping potential. In most applications 
of BEG, the harmonic potential is used nniin] 

1 f 7^a;2, d=l, 

(1.5) y(x) = -i + d=2, 

^ [ 7p^ + 7y+7zZ^ d=3, 

where jx, 7y and jz are three given positive constants. 

Define the energy functional 


( 1 . 6 ) E(4>) = [ 

jRd 




■ i/(x)|</>(x)|2 + |l<(>(x)|'^ - n4){-x.)L^4>{x.) 


dx. 


where / denotes the complex conjugate of /, then the ground state of a BEG is 
usually defined as the minimizer of the following nonconvex minimization problem 


(1.7) <j)g = arg min^gs E{(j)), 

where the spherical constraint S is defined as 


(1.8) 


S = 



< (X), 




It can be verified that the first-order optimality condition (or Euler-Lagrange equa¬ 
tion) of (1.7) is the nonlinear eigenvalue problem, i.e., find {p G M, 0(x)) such that 


= --V^</>(x) + y(x)</>(x) + /3|</>(x)|^(/)(x) - fJL^</>(x), X e 


(1.9) 

with the spherical constraint 

(1.10) ||<^f := / |<^(x)| 

jRd 


'dx = 1. 
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Any eigenvalue fi (or chemical potential in the physics literatures) of (L^-(l.lQ) can 
be computed from its corresponding eigenfunction 0(x) by [TOl [43l [45 


II = = E{(P) + f ^|^(x)|'^dx. 

jRd Z 


In fact, ( |1.9| ) can also be obtained from the GPE (1.3) by taking the anstaz '^(x, t) = 
0(x), and thus it is also called as time-independent GPE pill l4^ l45] . 

One of the two major problems in the theoretical study of BEG is to analyze and 
efficiently compute the ground state (pg in (1.7), which plays an important role in un¬ 


derstanding the theory of BEG as well as predicting and guiding experiments. Eor the 
existence and uniqueness as well as non-existence of the ground state under different 
parameter regimes, we refer to [T0l|39l|40l and references therein. Different numerical 
methods have been proposed for computing the ground state of BEG in the litera¬ 
tures, which can be classified into two different classes through different formulations 
and numerical techniques. The first class of numerical methods has been designed via 


the formulation of the nonlinear eigenvalue problem (1.9) under the constraint (1.10) 
with different numerical techniques, such as the Runge-Kutta type method 0133 ] for a 
BEG in ID and 2D/3D with radially/spherically symmetric external trap, the simple 
analytical type method [32] , the direct inversion in the iterated subspace method [48] , 
the finite element approximation via the Newton’s method for solving the nonlinear 
system m, the continuation method [25| and the Gauss-Seidel-type method [26j. In 


these numerical methods, the time-independent nonlinear eigenvalue problem (1.9) 


and the constraint (1.10) are discretized in space via different numerical methods. 


such as finite difference, spectral and finite element methods, and the ground state is 
computed numerically via different iterative techniques. The second class of numerical 
methods has been constructed via the formulation of the constrained minimization 
problem EH) with different gradient techniques for dealing with the minimization 
and/or projection techniques for handling the spherical constraint, such as the ex¬ 
plicit imaginary-time algorithm used in the physics literatures OlllSlEniEzlIlTj, the 
Sobolev gradient method [35] , the normalized gradient flow method via the backward 
Euler finite difference (BEED) or Eourier (or sine) pseudospectral (BEEP) discretiza¬ 
tion method 0 0 0 m m US] which has been extended to compute ground states 
of spin-1 BEG [l5l[T8], dipolar BEG [13] and spin-orbit coupled BEG [12], and the 
new Sobolev gradient method [30]. In these numerical methods, the time-independent 
infinitely dimensional constrained minimization problem EH is first re-formulated to 
a time-dependent gradient-type partial different equation (PDE) which is then dis¬ 
cretized in space and time via different discretization techniques and the ground state 
is obtained numerically as the steady state of the gradient-type PDE with a proper 
choice of initial data. 

Among those existing numerical methods for computing the ground state of BEG, 
most of them converge only linearly in the iteration and/or require to solve a large- 
scale linear system per iteration. Thus the computational cost is quite expensive 
especially for the large scale problems, such as the ground state of a BEG in 3D with 
an optical lattice potential or a rotating BEG with fast rotation and/or strong inter¬ 
action. On the other hand, over the last two decades, some advanced optimization 
methods have been developed for computing the minimizers of finite dimensional non- 
convex minimization problems, such as the Newton method via trust-region strategy 
[2a|42l[49] which converges quadratically or super-linearly. The main aim of this pa¬ 
per is to propose an efficient and accurate regularized Newton method for computing 
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the ground states of BEC by integrating proper PDE discretization techniques and 


advanced modern optimization methods. By discretizing the energy functional (1.6) 


and the spherical constraint (1.10) with either the finite difference, or sine or Eourier 


pseudospectral discretization schemes, we approximate the original infinite dimen¬ 


sional constrained minimization problem (1.7) by a finite dimensional minimization 
problem with a spherical constraint. Then we present an explicit feasible gradient 
type optimization method to construct an initial solution, which generates new trial 
points along the gradient on the unit ball so that the constraint is preserved auto¬ 
matically. The gradient type method is an explicit iterative scheme and the main 
costs arise from the assembling of the energy functional and its projected gradient 
on the manifold. Although this method often works well on well-posed problems, the 
convergence of the gradient type method is often slowed down when some parameters 
in the energy functional become large, e.g. /3 ^ 1 and O is near the fast rotation 
regime in (ini)- To accelerate the convergence of the iteration, we propose a regular¬ 
ized Newton type method by approximating the energy functional via its second-order 
Taylor expansion with a regularized term at each Newton iteration with the regular¬ 
ization parameter adjusted via the standard trust-region strategy [28l [42l [49]. The 
corresponding regularized Newton subproblem is a standard trust-region subproblem 
which can be solved efficiently by the gradient type method since it is not neces¬ 
sary to solve the subproblem to a high accuracy, especially, at the early stage of the 
algorithm when a good starting guess is not available. Eurthermore, the numerical 
performance of the gradient method can be improved by the state-of-the-art acceler¬ 
ation techniques such as Barzilai-Borwein steps and nonmonotone line search which 
guarantees global convergence [28l [42] 09] . In addition, we adopt a cascadic multigrid 
technique m to select a good starting guess at the finest mesh in the computation, 
which significantly reduces the computational cost. Extensive numerical experiments 
demonstrate that our approach can quickly reach the vicinity of an optimal solution 
and produce a moderately accurate approximation, even for the very challenging and 
difficult cases, such as computing the ground state of a BEC in 3D with an optical 
lattice potential or a rotating BEC with fast rotation and/or strong interaction. 

The rest of this paper is organized as follows. Different discretizations of the 
energy functional and the spherical constraint via the finite difference, sine and Eourier 
pseudospectral schemes are introduced in section In section we present the 
gradient type method and the regularized Newton algorithm for solving the discretized 
minimization problem with a spherical constraint. Numerical results are reported in 
section to illustrate the efficiency and accuracy of our algorithms. Einally, some 
concluding remarks are given in section Throughout this paper, we adopt the 
standard linear algebra notations. In addition, given X G the operators X, 

X*, 5R(X) and ^(X) denote the complex conjugate, the complex conjugate transpose, 
the real and imaginary parts of X, respectively. 

2. Discretization of the energy functional and constraint. In this section, 
we introduce different discretizations of the energy functional (1.6) and constraint 
(1.10) in the constrained minimization problem and reduce it to a finite dimen¬ 
sional minimization problem with a spherical constraint. Due to the external trapping 
potential, the ground state of (1.7) decays exponentially as |x| ^ oc [I0l|39l 00]. 


Thus we can truncate the energy functional and constraint from the whole space 
to a bounded computational domain U which is chosen large enough such that the 
truncation error is negligible with either homogeneous Dirichlet or periodic boundary 
conditions. We remark here that, from the analytical results uniisniiin], when 0 = 0, 
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i.e., a non-rotating BEC, the ground state <pg can be taken as a real non-negative 
function; and when 7^ 0, i.e., a rotating BEC, it is in general a complex-valued 
function, which will be adopted in our numerical computations. 


2.1. Finite difference discretization. Here we present discretizations of (1.6) 
and (1.10) truncated on a bounded computational domain U with homogeneous 
Dirichlet boundary condition by approximating spatial derivatives via the second- 
order finite difference (ED) method and the definite integrals via the composite trape¬ 
zoidal quadrature. Eor simplicity of notation, we only introduce the ED discretization 
in ID. Extensions to 2D and 3D without/with rotation are straightforward and the 
details are omitted here for brevity. 

Eor d = 1, we take U = (a, 6) as an interval in ID. Let h = {b — a)/N be 
the spatial mesh size with N a positive even integer and denote Xj = a ^ jh for 


j = 0,1 ,..., A^, and thus a = xq < xi <• - < xn-i < xn = b he the equidistant 
partition of U. Let <pj be the numerical approximation of for j = 0 , 1, ..., A/" 

satisfying 0o = 0(^o) = = 0(^Ar) = 0 and denote ^ = (^i,*** The 

energy functional ( |1.6[ ) with d = 1 and = 0 can be truncated and discretized as 


( 2 . 1 ) 


^(0) ^ / 

J a 

N-1 . 

i=o 

AT-l 


+ V{x)(j>{xf + ^4>{xf 


dx 


(x) + V(x)(t)(xf + ^(t){xY 


1 , 0i+l — + 0j_i 2 I ^/4 


dx 


2 N-l 


i=i ^ 

j=0 ^ '' j=l 


v{x^)4>] + 




= h 


JV-l 

i=i 




where A = (ajk) G i)x(Ar 1) is a symmetric tri-diagonal matrix with entries 

/ 

djk ^ 


7^ + V{xj), j — k, 

-2k’ \j-k\ = l, 

0, otherwise. 


Similarly, the constraint (1.10) with d = 1 can be truncated and discretized as 
(2.2) Uf ^ / 4>{xfdx = E / (t^{x)^dx ^h^4>]:= ||$||^ = 1, 

7=0 '^^3 7 = 1 


which immediately implies that the set S can be discretized as 

(2.3) Sh = {^e I Eh{^) < ||$||| = 1} . 

Hence, the original problem (EH with d = 1 can be approximated by the discretized 
minimization problem via the ED discretization: 


= arg 

5 


(2.4) 


















Denote Gh = VEh{^) be the gradient of Eh{^), notice ( |2.1[ ), we have 
(2.5) Gh := S/Eh{^) = 2h {A^ + , 

where G is defined component-wisely as for j = 1,..., A/" — 1. We 

remark here that, when the FD discretization is applied, the matrix ^4 is a symmetric 
positive definite sparse matrix. In addition, for the analysis of convergence and second 
order convergence rate of the above FD discretization, we refer the reader to [23l |52] . 


2.2. Sine pseudospectral discretization. For a non-rotating BEC, i.e. Cl = 
0, when high precision is required such as BEC with an optical lattice potential, 
we can replace the ED discretization by the sine pseudospectral (SP) method when 
homogeneous Dirichlet boundary conditions are applied. Again, we only present the 
discretization in ID, and extensions to 2D and 3D without rotation are straightforward 
and the details are omitted here for brevity. 

Eor d = 1, using similar notations as the ED scheme, similarly to (2.1), the energy 
functional (1.6) with d = 1 and D = 0 truncated on U can be discretized by the SP 
method as 


N-l r 


(2.6) 


E{cP) 


j = l 


-\4>j d^xx4>\x=xj 


where is the sine pseudospectral differential operator approximating the operator 
dxxi defined as 


(2.7) 


N-l 

^xx4^\x=Xj = ~ A 
1=1 


N 


= 1 , 2 , 


W-1, 


with the coefficients of the discrete sine transform (DST) of 4> G given 

as 


AT-l 


(2.8) = ( W ) ’ ^1= 


i=i 


b-a 


I = 1 , 2 ,--- ,N -1. 


Introduce V = diag(E(xi), • • • , V{xn-i)), A = diag(Ai,..., X%_i) and G = {cjk) G 
^(N-i)x(N-i) entries Cjk = sin for j,k = 1,..., — 1 and denote 4> = 

xT 

• • •, 0Ar-i j = Plugging (2.7) and (2.8) into (2.6), we get 


(2.9) 


E{(t)) ^ h 




AT-l 


i=i 


•=Ehm. 


where B G i)x(Ar i) is a symmetric positive definite matrix defined as 


( 2 . 10 ) 


B= —CAC + V. 


In fact, the first term in (2.9) can be computed efficiently at cost 0{N In N) through 
DST as 

N~ ~ 

(2.11) = T E + E 


1 = 1 


j = l 
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Again, the original problem (1.7) with d = 1 can be approximated by the discretized 
minimization problem via the SP discretization: 


( 2 . 12 ) 


$3 = arg 


Noticing (2.9), we have 


(2.13) Gh ;= VEhi^) = 2h (Bi> + ^i>3) =2h( + Vi> + 


2.3. Fourier pseudospectral discretization. For a rotating BEC, i.e. 7^ 0, 

due to the appearance of the angular momentum rotation, we usually truncate the 
energy functional ( |1.6[ ) and constraint (1.10) on a bounded computational domain U 
with periodic boundary conditions and approximate spatial derivatives via the Fourier 
pseudospectral (FP) method and the definite integrals via the composite trapezoidal 
quadrature. For simplicity of notation, we only introduce the FP discretization in 2D. 
Extensions to 3D are straightforward and the details are omitted here for brevity. 

For d = 2, we take U = [ai,6i] x [a2,&2] as a rectangle in 2D. Let hi = 
and ^2 = be the spatial mesh sizes with and N 2 two positive integers 

and denote Xj = ai ^ jhi for j = 0, 1 , ..., A^i, yk = ^2 + kh 2 for k = 0, 1 , ..., A^2- 
Denote h = max{hi, ^2} and Ujk = {xj^Xjj^i) x {yk^yk+i)- Let be the numerical 
approximation of (j){xj^yk) for j = 0, 1 ,..., A^i and /c = 0, 1 ,..., A^2 satisfying (t)jN 2 = 
(j)jn for j = 0,1 ,..., and (j)N^ k = ^o/c for /c = 0,1 ,..., A^2 and denote ^ G 

C^ix^2 energy functional (1.6) with d = 2 can be truncated and discretized as 


E{ctp) 


ohi ph 2 r 


-\4>^4> + V{x,y)\4>\'^ + 


■ iflpixdy -ydx)4> 


dxdy 


n 

J a\ Ja 
N\ — 1 N 2 — 1 p 

j=0 /c=0 ■- 

a ^2 r _ /i 1 

^1^2 F F ( 2 dL^ljk + 2 ^yy^\jk + 

j=0k=0 *- ^ 


-^(f>A(p + V{x,y)\(f>f + ^\4>\^ + in4>{xdy - ydx)(p 


dxdy 


(2.14) 

+V{xj,yk)\<pjk\^ + 

^\4>jk\^ ajk ■■= 

-- Eh{^), 

where 

f 1 

l<3<N^- 

1, l<k<N2-l, 


ajk = s 1/4 

j = O&A: = 0, N 2 or j = Ni^k = 0, A' 2 , 


1 1/2 

otherwise. 


and the Fourier pseudospectral differential operators are given as 


Ni/2-1 


Ni/2-1 




dL^\jk = - F 

(2.15) 

p=-iVi/2 

p=-iVi/2 

N 2 I 2 -I 

N 2 I 2 -I 


8,''^U= E 

. 7(2) 

^Vq<Pjqe , 

9ly^\jk = - E 


9=-JV2/2 

q=-N2/2 
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with 


(2.16) 


1 

^pk ~ j\j- 

j=0 




N 2-1 




k=0 


27 r jp 
. * N-, 

Xp — 

2 'kp 

P = 

iVi 

iVi 

-1, 


bi-aP 

2 

2 

• 2TTkq 
: «2 , 

Vq = 

27 rq 


N2 

N2 

- 1. 

62 — ^2 ’ 

Q. = 

2 

■’ 2 


Plugging (2.15) and (2.16) into (2.14), the discretized energy functional Eh{^) can 
be computed efficiently via the fast Fourier transform (FFT) as 


Eh{^) = hih2 


N2 iVi/2-1 / .2 \ 

gai.ATi 

k=0 p=-Ni/2 \ / 


Ni N 2 / 2 -I / 2 

Y^ajlN2 I ^ - XjTJqil ) \4>'. 

j =0 q=-N 2 l 2 

Ni N 2 


( 2 ) 12 
IQ I 


(2.17) 


-\-hlh2 EE O^jk 

j=0 k=0 


2 I |4 


^{^j^yk)\4>jk\ + r^l^jkl 


Similarly, the constraint (1.10) with d = 2 can be truncated and discretized as 

Nx — 1 N 2 — 1 


( 2 . 18 ) Uf ^ f f \(t){x,y)\'^dxdy hih 2 Y E 

dai Ja2 

which immediately implies that the set S can be discretized as 

(2.19) Sh = {^e I Eh{^) < 00 , ||$||2 = 1} . 

Hence, the original problem with d = 2 can be approximated by the discretized 
minimization problem via the FP discretization: 


= 11 ^ 11 ^ = 1 , 


( 2 . 20 ) 


= arg £;,»($). 


Noticing (2.17), similarly to (2.13), Gh = \/Eh{^) can be computed efficiently via 
FFT in a similar manner with the details omitted here for brevity. 

3. A regularized Newton method by trust-region type techniques. It is 

easy to see that the constrained minimization problems (2.4), (2.12) and (2.20) can 
be written in a unified way via a proper rescaling 


(3.1) 


1 ^ 

Xg := arg := -XMX + , 


i=i 


where M is a positive integer, o is a given real constant, A G is a Hermitian 

matrix and the spherical constraint is given as 


Em 


|x = (Xi,X2,...,Xm)^ ec^ I ||x||i = i| 
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We first derive the optimality conditions of the problem (3.1). The gradient and 
Hessian of J-{X) can be written explicitly. 

Lemma 3.1. The first and second-order directional derivatives of T{X) along a 
direction D G are: 


M 


(3.2) 


VT{X)[D] = ^{D*AX) + ia'^iXjXV^XjDj), 


i=i 


M 


(3.3) X^X{X)[D,D]=D*AD + 4aJ2 [(X^XViD^Dj)+2?ft{XjDV^^ 

i=i 


Define the Lagrangian function of (3.1) as 
(3.4) 


L{X,d)=X{X)--{\\Xr2-l), 


then the first-order optimality conditions of 


are 


(3.5) 

(3.6) 


G-8X = 0, 

IAI|2 = i, 


where G = VX{X) is the gradient of X{X). Multiplying both sides of (3.5) by X* 
and using (|3.6|), we have 0 = X*G. Therefore, (|3.5|) becomes 


(3.7) {I-XX^)G = A{X)X = 0, with M(X) = GX* - XG*. 

By definition, A{X) is skew-symmetric at every X. 

By differentiating both sides of X*X = 1, we obtain the tangent vector set of the 
constraints: 


(3.8) 


Tx--{Z : X*Z = 0}. 


The second-order optimality conditions is described as follows. 

Lemma 3.2. 1) (Second-order necessary conditions, Theorem 12.5 in M) Sup¬ 
pose that X G is a local minimizer of the problem (3.1). Then X satisfies 


(3.9) V^X(X)[D,D] - OD^D >0, VD G where 0 = VX(X)*X. 


2) (Second-order sufficient conditions, Theorem 12.6 in M) Suppose that for 
X G , there exists a Lagrange multiplier 0 such that the first-order conditions are 
satisfied. Suppose also that 


(3.10) 


V‘^X{X)[D,D]-0D*D > 0, 


for any vector D G Tx- Then X is a strict local minimizer for (3.1). 


3.1. Construct initial solutions using feasible gradient type methods. 

In this subsection, we consider to solve the problem by following the feasible 

method proposed in m- The description of the algorithm is included to keep the 
exposition as self-contained as possible. Observe that A{X)X is the gradient of X(X) 
at X projected to the tangent space of the constraints. The steepest descent path 
is Y{r) := X — tA{X)X, where r is a positive constant representing the step size. 
However, this Y(r) does not generally have a unit norm. 
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An alternative implicit updating path is 


(3.11) Y{t) —X- tA{X){X + y(T)) ^ Y{t) = (/ + TA{X)y^ {/ - tA{X)) X. 


Then the fact that (/ rA{X))~^ {I — rA{X)) is orthogonal for any r > 0 gives 

||h ^('?")||2 = ll-^lb = 1, he., the constraints are preserved at every r. The closed-form 
solution of Y (r) can be computed explicitly as a linear combination of X and G, in 
which the linear coefficients are determined by r, ||-^||2 5 ||G ||2 and X*G. 

Theorem 3.3. For every r > 0, Y{r) of ( 3.11| ) satisfies ||T(r )||2 = hn 

addition, Y{r) is given in the elosed-form as 


(3.12) 


Y{r) = a{r)XYp{r)G, 


where 

^ {l+TX*Gf-y\\X\\l\\G\\l ^ -2 t\\X\\1 

i-Tyx*Gf + Tyx\\i\\G\\V > i-Tyx*G)^ + Tyxg\\G\\r 


We refer to m for the details of the proof of this theorem. 

A suitable step size r can be chosen by using a nonmonotone curvilinear (as our 
search path is on the manifold rather than a straight line) search with an initial step 
size determined by the Barzilai-Borwein (BB) formula [20]. They were developed 
originally for the vector case in [20]. At iteration k, the step size is computed as 


(3.13) 




tr ((5(fe-i))*5(fe-i)) 

|tr((5(*^-i))*iy(*-i))| 




|tr ((5(fe-i))*vi/('=-i)) I 
tr((lF(''-i))*lF(''-i)) 


where and = A{X^^^)X^^^ -A{X^^-^^)X^^-^\ When 

or is not bounded, they are reset to a finite number. 

In order to guarantee convergence, the final value for is a fraction of or 


^/c,l 


determined by a nonmonotone search condition. Let T(r) be defined by (3.11) 
^( 0 ) _ x(X*^^^), + 1 and = 1. The new points are generated 

iteratively in the form := T(^)(r^^^) with 

Here m is the smallest nonnegative integer satisfying 


r(^) = or 


(3.14) 


X(T(^)(t(^))) < G^^^ -piT^^^||^(X^^))X^^)||2, 


where each reference value is taken to be the convex combination of G^^^ and 

X(X^^+^^) as F In Algorith m [l] below, we 

specify our method for solving the constrained minimization problem (3.1) obtained 
from the discretization of the ground state of BEG. Although several backtracking 
steps may be needed to update the X^^+^\ we observe that the BB step size or 
is often sufficient for ( |3.14[ ) to hold in most of our numerical experiments. 

We can establish the convergence of Algorithm as follows. 

Theorem 3.4. Let {X^^^ \ k > Y\ he an infinite sequenee generated by the 
Algorithm^ Then either ||^(X^^^)X^^^ ||2 = 0 for some finite k or 


liminf ||^(X^^^)X ^^^||2 = 0. 

k^oo 
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Algorithm 1: A feasible gradient method 


1 

2 

3 

4 

5 

6 


Given set pi,G (0,1), /c = 0. 
while stopping conditions are not met do 

Compute ^ or ^ where m is the smallest 

nonnegative integer satisfying the condition (|3.14|). 

Set ^ Y{r). 

Q(/c+i) ^ ^Qik) ^ I ^(/c+i) ^ 

k ^ k^l. 


Proof. Since the energy function P{X) is differentiable and its gradient VP{X) is 
Lipschitiz continuous, the results can be obtained using the proofs of in a similar 
fashion. □ 

Remark 3.5. The convergence of the full sequence can be ensured if a 

monotone line search is used. Given d > 0,pi,(5 G (0,1), the Armijo point at X^^^ 
is defined as where Y{t) is the curve (|^ TT] ), and m is the 

smallest nonnegative integer satisfying 


(3.15) 




Using the proofs of Theorem j.S.l and Corollary 4.3.2 m in a similar fashion, we 
can prove that lim/c^oo \\A{X^^^)X ^^^\2 = 0. 

3.2. A regularized Newton method for computing ground states of 

BEC. In general, the Algorithm works well in the case of weak interaction and 
slow rotation, i.e. \I3\ and \ft\ are small in the energy functional ( |1.6| ). However, its 
convergence is often slowed down in the case of strong interaction and/or fast rota¬ 
tion, i.e., when one of the parameters becomes larger, and thus it can take a lot of 
iterations to obtain a highly accurate solution. Usually, fast local convergence can¬ 
not be expected if only the gradient information is used, in particular, for difficult 


non-quadratic problems. Observe that the most difficult term in (3.1) is the quartic 
function \Xi\^. A Newton method is to replace J-{X) by its second-order Taylor ex¬ 
pansion. In order to ensure the global convergence of the Newton’s method, we adopt 
the trust region methc 
surrogate function as: 


the trust region method [SHiiiaiin] by adding a proximal term ||A — in the 


1 

W^^\X) := VP{X^^^)[X-X^^^]P-V‘^P{X^^^)[X-X^^\X-X^^^]P—\\X-X^^^\\1, 


where 5^^^ > 0 is a regularization parameter. Using Lemma [3^ we obtain that 

W^^\x) = W^^\x) + constant. 


where 


= -X*AX + 4a 


N 


+ 2a 


N r 

E ( 

j = l ^ 


j = l 




\Xj - + 25R (y {Xj - xf'^ 


2 
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The gradient of W^^\X) is 


(VVFW(X))^- = 

We next present the regularized Newton framework starting from a feasible initial 
point and the regularization parameter . At the k-th iteration, our regularized 
Newton subproblem is defined as 


(3.16) 


mm 

IAII2=i 


lyW(X) 


The subproblem (3.16) is the so-called trust-region subproblem. Since the dimension 
M in is usually very large so that the discretization error of ( |1.7[ ) can be small, 
the standard algorithms for solving the trust-region subproblem [28l HU 09] usually 
cannot be applied to (3.16) directly. Hence, we still use a gradient-type method 
similar to the one described in subsection 3.1 to solve (3.16). The method is ideal for 


solving these regularized Newton subproblems since it is not necessary to solve these 
subproblems to a high accuracy, especially, at the early stage of the algorithm when 
a good starting guess is not available. 

Let be an optimal solution of ( |3.16[ ). Generally speaking, an algorithm 

cannot be guaranteed to converge globally if is set directly to the trial point 

obtained from a model with a fixed In order to decide whether the trial 

point should be accepted and whether the regularization parameter should be 
updated or not, we calculate the ratio between the actual reduction of the objective 
function X{X) and predicted reduction: 


(3.17) 






If > 0, then the iteration is successful and we set = Z^^^; otherwise, 

the iteration is not successful and we set = X^^\ that is, 


(3.18) 


= 


Z^^\ if > 771 , 
X^^\ otherwise. 


Then the regularization parameter is updated as 

( (0,(5^^^], if p^^^ > 7 / 2 , 

(3.19) ^(^+1) ^ J if Pi < p^^^ < 772 , 

[ 72(5^^^], otherwise. 

where 0 < 771 < 772 < 1 and 1 <71 <72- These parameters determine how aggres¬ 
sively the regularization parameter is decreased when an iteration is successful or it 
is increased when an iteration is unsuccessful. In practice, the performance of the 
regularized Newton algorithm is not very sensitive to the values of the parameters. 

The complete regularized Newton algorithm to solve is summarized in the 
Algorithm 

The convergence of the Algorithm can also be established as follows. 

Theorem 3.6. Let {X^^^ : k > 0} be an infinite sequenee generated by the 
Algorithm^ Then either ||^(X^^^)X^^^ II 2 = 0 for some finite k or 

lim ||^(X^^^)X ^^^||2 = 0. 

k^oo 
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Algorithm 2: A regularized Newton method 


1 

2 

3 

4 

5 

6 

7 

8 


Given a feasible initial solution with ||X ^^^||2 = 1 and initial 
regularization parameter >0. Choose 0 < < 772 < 1, 1 < 7i < 72 - 

Call Algorithm to minimize problem (3.1) to a certain low accuracy for a 
feasible solution Set iteration k := 1. 

while stopping conditions are not met do 

Solve (3.16) to obtain a new trial point . 

Compute the ratio via (3.17). 

Update from the trial point based on (3.18). 

Update 6^^^ according to (3.19). 
k ^ k ^1. 


Proof. Since the energy function J-{X) is differentiable and its gradient XJ-{X) is 
Lipschitiz continuous, the results can be obtained using the proofs of [50] in a similar 
fashion. □ 

The discretization of (EH) on a fine mesh usually leads to a problem of huge size 
(M ^ 1) whose computation cost is very expensive, especially for high dimensional 
case. A useful technique is to adopt the cascadic multigrid method [21], i.e. solve 
the minimization problem (1.7) on the coarsest mesh, and then use the obtained 
solution as the initial guess of the problem on a fine mesh, and repeat until we obtain 
the solution on the finest mesh. We present the mesh refinement technique via the 
cascadic multigrid method in the Algorithm where the discretized problems are 
solved from the coarsest mesh to the finest mesh. 


Algorithm 3: A cascadic multigrid method for mesh refinement 

1 Given an initial mesh and set k = 0. 

2 while convergence is not met do 

3 


Use as an initial guess on the kth mesh to calculate the optimal 
solution of the minimization problem (3.1) using the Algorithm]^ 

Refine the mesh uniformly to obtain 
k ^ k 1. 


4. Numerical results. In this section, we report several numerical examples to 
illustrate the efficiency and accuracy of our method. All experiments were performed 
on a PC with a 2.3GHz CPU (i7 Core) and the algorithms were implemented in 
MATLAB (Release 8.1.0). In our experiments, the Algorithmis called to compute 
the ground state of non-rotating BEC, i.e., U = 0, since it is a relatively easy problem. 
The algorithm is stopped either when a maximal number of K iterations is reached 
or when 


(4.1) 


||X(fe+i) -XW||oo 

q-(k) 


< ^ 0 - 


The default values of Cq and K are set to be 10 ^ and 2000, respectively. In order to 
test the spectral accuracy of the SP discretization, a tighter stopping criterion is taken. 
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A normalization step is executed if |X*X — 1| > to enforce the feasibility. For 

non-rotating BEC with strong interaction, i.e., /3 ^ 1, the initial solution is usually 
chosen as the Thomas-Fermi (TF) approximation [T0l|T6l[45] 


(4.2) 


00 (x) = 


— V (x) 


where | (^) ’ ^ d = 1, 2 and 3, respectively. 

Since the Algorithm may converge slowly for computing the ground state of rotating 


.0. 

1/2 


if F(x) 
otherwise, 


BEC, i.e., O 7 ^ 0, we choose the regularized Newton method (i.e.. Algorithm 
together with the cascadic multigrid method for mesh refinement (i.e.. Algorithm 
and it is terminated when 


3) 


(4.3) 


||X(''+l) -X('')||oo <<5o, 


where the default value of (5o is set to 10“^. Let (pg be the “exact” ground state 
obtained numerically with a very fine mesh and we denote its energy and chemical 
potential as Eg = E{(j)g) and jig = /i(0^), respectively. To quantify the ground state, 
one important quantity is the root mean square which is defined as 


(4.4) 


CCr 


= \\a(t>g\\mu) = ^ 


a = x^y or z. 


4.1. Accuracy test and results in ID. We take d = 1 and = 0 in (1.7) and 


(1.6) and consider two different trapping potentials 

Case I. A harmonic oscillator potential ( |1.5[ ) with d = 1, = 1 and [3 = 400. 

Case II. An optical lattice potential V{x) = ^ E 25sin^(i^) and 0 = 250. 

The ground state is numerically computed by the Algorithm on a bounded 
computational domain U = (—16,16) which is partitioned equally with the mesh 
size h. In ord er to compare the accuracy of the FD and SP discretizations, we set 


eo = 10 in (4.1). Let 0™ and be the numerical ground states obtained with 


^9,h 


the mesh size fi by using FD and SP discretization, respectively. Table depicts the 
numerical errors for Case I, and respectively. Table [^for Case 11. 


Table 1 

Accuracy of the FD and SP discretizations for Case I in 


Mesh size 

h = l 

h = 1/2 

h = 1/4 

h = 1/8 

max \4)g - 

\Eg-E{&\ 

K-MOl 

2.06E-03 

8.59E-04 

2.21E-02 

1.24E-03 

2.66E-04 

9.48E-05 

2.88E-04 

6.46E-05 

3.49E-05 

7.43E-05 

1.59E-05 

8.60E-06 

max \4)g - 
\Eg-E{<l>p\ 

1.31E-03 

5.69E-05 

1 .66E-02 

7.04E-05 

2.64E-06 

8.71E-05 

1.95E-08 

8.45E-12 

9.55E-10 

5.01E-13 

2.17E-13 

2.52E-12 


From Tables and it is observed that the SP discretization is spectrally accu¬ 
rate, while the FD discretization has only second order accuracy for computing the 
ground state of BEC in ID. Hence, when high accuracy is required, the SP discretiza¬ 
tion is preferred since it needs much fewer grid points, and thus it saves significantly 
memory cost and computational cost. 
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Table 2 

Accuracy of the FD and SP discretizations for Case II in 


Mesh size 

h = l 

h = 1/2 

h = 1/4 

h = 1/8 

max \4,g - (^>71 

\F-EACh)\ 

K-MOl 

1.02E-02 

2.66E-02 

1.27E-01 

5.81E-03 

8.39E-03 

4.05E-03 

9.97E-04 

2.03E-03 

8.28E-04 

2.50E-04 

5.02E-04 

2.08E-04 

max 

7.98E-03 

4.22E-04 

9.76E-02 

1.21E-03 

1.96E-04 

4.11E-03 

2.22E-06 

4.99E-08 

5.61E-07 

1.90E-11 

7.53E-13 

9.17E-13 


For comparison with existing numerical results in the literatures [annnHiiiin], 
Figurej^plots the ground states obtained by the SP discretization for cases I and 

IL In addition, their energy, chemical potential and root mean squares are obtained 
as for Case I: Eg = 21.3601, jig = 35.5775 and Xj-ms = 3.7751; and for Case II: 
Eg = 26.0839, IIg = 38.0692 and Xrms = 3.3609. These numerical results agree very 
well with those reported in the literatures m uni m USE]. 

Fig. 1. Ground states (fgix) for Case I (left) and Case II (right) in 




4.2. Accuracy test and results in 3D. We take d = 3 and O = 0 in (1.7) and 


(1.6) and consider two different trapping potentials m- 

Case I. A harmonic oscillator potential (1.5) with d = 3, = 1, = 2, 7 ;^ = 4 

and p = 200. 

Case II. A harmonic oscillator potential and a potential of a stirrer corresponding 
to a far-blue detuned Gaussian laser beam 

V{x, y, z) = l{x^+ ly + Ay) + 

with y^ = 1, 7;^ = 2, cjo = 4, d = 1, ro = 1 and P = 200. 

Again, the ground state is numerically computed by the Algorithm on bounded 
computational domains U = (—8,8) x (—6,6) x (—4,4) and U = (—8,8)^ for Case I 
and II, respectively, which are partitioned uniformly with the same number of nodes 
in each direction. Let h be the mesh size in the x-direction. Again, we set cq = 10“^^ 
in (4.1). Let 0™ and be the numerical ground states obtained with the mesh 
size ti by using FD and SP discretization, respectively. Table depicts the numerical 
errors for Case I, and respectively. Table [^for Case IL 
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Table 3 

Accuracy of the FD and SP discretizations for Case I in 


Mesh size 

II 

to 

h = l 

h = 1/2 

h = 1/4 

max \4)g - 

K-MOl 

2.28E-02 

1.26E-01 

4.45E-02 

5.16E-03 

5.82E-02 

3.10E-02 

l.llE-03 

1.44E-02 

9.40E-03 

2.51E-04 

3.41E-03 

2.23E-03 

max 

K-MC)I 

l.lOE-02 

l.OlE-01 

1.57E-01 

1.68E-03 

6.49E-05 

4.17E-03 

8.68E-06 

1.45E-08 

5.48E-07 

7.34E-10 

1.09E-11 

1.55E-11 

Table 4 

Aeeuraey of the FD and SP diseretizations for Case II in 

U-2. 

Mesh size 

h = 2 

h = 1 

h = 1/2 

h = 1/4 

max \(Pg - 
\Eg-E{4>Jf^)\ 

K-MOl 

1.61E-02 

6.76E-01 

5.37E-01 

7.92E-03 

6.06E-02 

6.16E-02 

1.69E-03 

1.33E-02 

8.09E-03 

3.92E-04 

3.16E-03 

1.92E-03 

max 103 - 

\E,-E{4g^)\ 

1.69E-01 

1.87E-01 

5.69E-01 

2.57E-03 

6.69E-03 

2.21E-02 

4.38E-05 

9.55E-06 

7.79E-06 

1.18E-08 

6.34E-12 

9.85E-11 


Again, from Tables and it is observed that the SP discretization is spectrally 
accurate, while the FD discretization has only second order accuracy for computing 
the ground state of BEC in 3D. Hence, when high accuracy is required and/or the 
solution has multiscale phenomena, the SP discretization is preferred since it needs 
much fewer grid points, and thus it saves significantly memory cost and computational 
cost. 

Again, for comparison with existing numerical results in the literatures 
[161 Hr], Figure]^ plots the ground states 0 ^(t,O, 2 ;) obtained by the SP discretization 
for cases I and II. In addition, their energy, chemical potential and root mean squares 
are obtained as for Case I: Eg = 8.3345, jig = 11 . 0102 , Xj-ms = 1.6710, ^rms = 0.8751, 
and Zrms = 0.4884; and for Case II: Eg = 5.2696, = 6.7019, x^ms = 1-3744, 

^rms = 1.4358 and Zrms = 0.7043. These numerical results agree very well with those 
reported in the literatures laiioiiiiiiiiiiT]. 

To demonstrate the high resolution of the SP discretization and compare our 
algorithm with existing numerical methods Huniiiii, we also apply our algorithm 
to compute the ground state of BEC in 3D with a combined harmonic and optical 
lattice potential m as 


(4.5) V{x,y,z) 





+ sin2 



together with different interaction constants (3 = 100, 800 and 6400. The ground state 
is numerically computed by the Algorithm on bounded computational domains 
U = (— 8 , 8 )^ for P = 100 and 800, and U = (— 12 , 12 )^ for /3 = 6400, which are 
partitioned uniformly with the number of nodes = A ^2 = = 2 ^ + 1 in each 

direction. The stopping criterion is set to the default value. 

Table depicts the maximum value of the wave function max|0^p, the energy 
E{(j)g)^ the chemical potential jig and the root mean squares Xmis, ^rms and Zrms for 
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Fig. 2. Ground states 0^(a:,O, 2 ;) for Case I (left) and Case II (right) in ^4-^- 



Table 5 

Comparison of numerieal results eomputed by Algorithm^ (top half part with rows 2-4 ) and 
BESP implemented in GPELah (bottom half part with rows 5-1) for trapping potential with 

different /3. 


/3 

max 10^1^ 

E{ci>g) 

ftg 

Trms 

^rms 

^rms 

iter 

nfe 

cpu (s) 

100 

0.2536 

23.2356 

27.4757 

1.8716 

1.8716 

1.8716 

112 

115 

76.47 

800 

0.0490 

33.8023 

40.4476 

2.6620 

2.6620 

2.6620 

260 

279 

183.34 

6400 

0.0098 

52.4955 

63.7146 

3.3685 

3.3685 

3.3685 

305 

327 

217.03 

100 

0.2536 

23.2356 

27.4757 

1.8717 

1.8717 

1.8717 

188 

- 

914.53 

800 

0.0490 

33.8023 

40.4476 

2.6620 

2.6620 

2.6620 

494 

- 

2513.75 

6400 

0.0098 

52.4955 

63.7149 

3.3684 

3.3684 

3.3684 

747 

- 

4014.17 


different interaction constants /d, as well as the number of iterations (iter), the number 
of function evaluations (nfe) and the computational time (cpu). For comparison, we 
also display numerical results obtained by using the BESP method implemented in 
GPELab [7] (a MATLAB toolbox designed for computing ground state and dynamics 
of BEG) with time step taken as At = 10“^ and all other setting the same as above. 
In addition, Eigurej^ shows the isosurface plots and their corresponding slice views of 
the ground states for different /3. 

Erom Table we can see that the Algorithm converges to the ground state 
much faster than the BESP method presented in [TUllTi] for all p in computing the 
ground state of BEG in 3D. 


4.3. Results for rotating BEC in 2D. We take d = 2 and the harmonic 
potential with 7 a, = = 1 in ( |1.7[ ) and ( |1.6[ ) and consider different [3 and 

D. The ground state is numerically computed by the regularized Newton method 
(i.e. Algorithm with the EP discretization on bounded computational domains 
U = (—10,10)^ and U = (—12,12)^ for p = 500 and P = 1000, respectively. The 
domains are partitioned uniformly with the number of nodes Ni = N 2 = 2^ -\- 1 in 
each direction. In our computations, in the Algorithm we first call the gradient 
type method, i.e.. Algorithmwith a maximum number of iterations i^init = 100 to 
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Fig. 3. Isosurfaces (left column) and their eorresponding sliee views (right column) o f the 
ground states (l)g{x^y,z) of a BEC in 3D with eombined harmonie and optieal lattiee potential 
for different j3 = 100, 800, 6400 (from top to bottom). 




obtain a good initial guess Then the regularized Newton subproblem is solved 

by the Algorithm up to a maximum number of iterations K^uh = 200. In order to 
reduce the computational cost, the cascadic multigrid method (i.e., Algorithm]^ is 
applied for mesh refinement with the coarsest mesh chosen with the number of 
nodes A^i = A ^2 = 2^ + 1 in each direction. 

For a rotating BEC, the ground state is a complex-valued function, and thus 
it is very tricky to choose a proper initial data such that the numerical result is 
guaranteed to be the ground state. Similarly to those in the literatures m, here we 
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test our algorithms with the following different initial solutions 


(a) (t)a{x,y) = 


J_e-(®^+y^)/2 


(b) (f>b{x,y) 


x + iy 

A 


(6) (pi{x,y) = (pb{x,y), 


(c) Mx,y) = 

(d) (t)d{x,y) = 


[4>a{x,y)) + (f>b{x,y)\/2 
\\[4>a{x,y)) + 4>b{x,y)]/2\y 
(1 - fl)(pa{x, y)) + 9.4>b{x, y) 
11(1 - Q.)4)a{x,y)) + Q.4)b{x,y)\\ 


(c) (t)c{x,y) = 4>c{x,y), 

, {d) (t)b{x,y) = ^d{x,y). 


Table displays the energy obtained numerically with different initial data se- 
lected in the above with p = 500 for different Q. = 0.00, 0.25, 0.50, 0.60, 0.70, 0.80, 
0.90 and 0.95 (in the table, we use a “f” sign to indicate the one with the lowest 
energy among different initial data for given [3 and f2), and Table summarizes the 
lowest energy among different initial data and the corresponding number of iterations 
and computation time for P = 500 with different ft. Figure plots the ground state 
density for P = 500 with different ft. In addition. Tables [8p| and Figure]^ 

present similar numerical results for p = 1000. 


Table 6 

Energy obtained numerieally with different initial data of rotating BECs for (3 = 500 and 
different Cl in ^4-3- 


a 

0.00 

0.25 

0.50 

0.60 

0.70 

0.80 

0.90 

0.95 

(a) 

8.5118 

8.5118 

8.0246 

7.5890 

6.9731 

6.1016 

4.7778 

3.7417 

(b) 

8.5118 

8.5106 

8.0246 

7.5845 

6.9731 

6.1055 

4.7778 

3.7417 

(b) 

8.5118 

8.5118 

8.0197+ 

7.5890 

6.9731 

6.1016 

4.7778 

3.7416 

(c) 

8.5118 

8.5106 

8.0246 

7.5890 

6.9726 

6.1016 

4.7778 

3.7417 

(c) 

8.5118 

8.5118 

8.0246 

7.5890 

6.9731 

6.0997 

4.7778 

3.7415 

(d) 

8.5118t 

8.5106t 

8.0246 

7.5890 

6.9726t 

6.0997+ 

4.7778t 

3.7415t 

(d) 

8.5118 

8.5118 

8.0246 

7.5845+ 

6.9731 

6.1016 

4.7778 

3.7416 


Table 7 

Ground state energy, the number of iterations for the regularized Newton method (iter) on the 
finest mesh and the total eomputational time (epu) of rotating BECs for /3 = 500 and different Cl in 

U.3. 


n 

0.00 

0.25 

0.50 

0.60 

0.70 

0.80 

0.90 

0.95 

iter 

3 

3 

3 

128 

49 

18 

69 

4 

epu (s) 

1.14 

18.71 

41.57 

355.63 

147.03 

130.87 

286.12 

56.08 

Energy 

8.5118 

8.5106 

8.0197 

7.5845 

6.9726 

6.0997 

4.7778 

3.7415 


From Tables among those different initial data, either (d) or (d) gives the 
lowest energy in most cases. Thus, in practical computations, we recommend to choose 
either (d) or (d) as the initial data. Also, it is observed that the regularized Newton 
algorithm converges quickly to the stationary solution within very few iterations, even 
for strong interaction, i.e., P ^ and fast rotation i.e., ft is near 1. Compared with 
the normalized gradient flow method via BEFD or BESP discretization for computing 
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Fig. 4. Plots of the ground state density |0^(a:,y)p - eorresponding to the energy listed in the 
Table^ - of rotating BECs for j3 = 500 and different tl in %4--3- 
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Table 8 

Energy obtained numerically with different initial data of rotating BECs for (3 = 1000 and 
different Cl in §4-3- 


£2 

0.00 

0.25 

0.50 

0.60 

0.70 

0.80 

0.90 

0.95 

(a) 

11.9718 

11.9718 

11.0954t 

10.4392 

9.5335 

8.2610 

6.3608 

4.8830 

(b) 

11.9718 

11.9266 

11.1326 

10.4392 

9.5283 

8.2610 

6.3607 

4.8825 

(b) 

11.9718 

11.9266 

11.1054 

10.4392 

9.5335 

8.2631 

6.3607 

4.8827 

(c) 

11.9718 

11.9165 

11.1054 

10.4392 

9.5289 

8.2610 

6.3607 

4.8823^ 

(c) 

11.9718 

11.9165 

11.1326 

10.4392 

9.5283 

8.2610 

6.3607 

4.8825 

(d) 

11.9718 

11.9266 

11.1054 

10.4392 

9.5289 

8.2632 

6.3608 

4.8825 

(d) 

11.9718t 

11.9165^ 

11.1326 

10.4392t 

9.5283t 

8.2610t 

6.3607t 

4.8830 
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Table 9 

Ground state energy, the number of iterations for the regularized Newton method (iter) on the 
finest mesh and the total eomputational time (cpu) of rotating BECs for (3 = 1000 and different Q. 
in ^4-3- 


n 

0.00 

0.25 

0.50 

0.60 

0.70 

0.80 

0.90 

0.95 

iter 

3 

3 

3 

10 

10 

72 

41 

157 

cpu (s) 

1.18 

28.52 

108.98 

106.86 

105.28 

313.67 

825.12 

751.72 

Energy 

11.9718 

11.9165 

11.0954 

10.4392 

9.5283 

8.2610 

6.3607 

4.8823 


Fig. 5. Plots of the ground state density \(f)g(x,y)\^ - eorresponding to the energy listed in the 
- of rotating BECs for (3 = 1000 and different Q in %4-3- 

Q=0.25 
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ground state of a rotating BEC 13 la [a SOI ss], the regularized Newton algorithm 
significantly reduces the computational time. 


4.4. Application to compute asymmetric excited states. When the trap¬ 
ping potential E(x) in (1.7) is symmetric and the BEC is non-rotating, similarly to 
those numerical methods presented in the literatures [loi [HI [13 na, our numerical 
methods can also be applied to compute the asymmetric excited states provided that 
the initial data is chosen as an asymmetric function. To demonstrate this, we take 
d = 2, = 0 and P = 500 in ( |1.7[ ) and the trapping potential is chosen as a combined 

harmonic and optical lattice potential 


(4.6) 


y{x,y) = i +y^) +50 


sm' 


If) 


■ sm 


If) 


The ground and asymmetric states are numerically computed by the Algorithm via 
the SP discretization on the bounded computational domain U = (—16,16)^ which is 
partitioned uniformly with the number of nodes = A 2 = 2^ + 1 in each direction. 


The initial data is chosen as the TE approximation (4.2) for computing the ground 


state <pg^ as for the asymmetric excited state in the x- 

direction 0io, as (j){){x^y) = for the asymmetric excited state in the 

^-direction and as (^^{x^y) = for the asymmetric excited state in 

both X- and d-birections ^n, respectively. The stopping criterion is set to the default 


value. Table [TQ] lists different quantities of these states and computational cost by our 
algorithm. In addition, Eigurej^ shows contour plots of these states. 

Table 10 

Different quantities of the ground and asym metri c excited states and the corresponding compu¬ 
tational cost for a BEC in 2D with the potential Hi and [3 = 500 in §4-4- 


0 

max 0 


m(<(') 

Trms 

^rms 

iter 

nfe 

cpu (s) 

<Pg 

0.0820 

32.2079 

41.7854 

2.9851 

2.9851 

365 

380 

3.99 

010 

0.0746 

34.6053 

43.8248 

3.3029 

2.8741 

285 

301 

3.18 

001 

0.3749 

34.6053 

43.8248 

2.8741 

3.3029 

272 

288 

3.03 

011 

0.0666 

37.0864 

46.1442 

3.1434 

3.1434 

117 

125 

1.32 


Erom Table IT and Eigure we can see that our algorithm can be used to 
compute the asymmetric excited states provided that the initial data is taken as 
asymmetric functions. The numerical results from our algorithm agree very well with 
those reported in the literatures [13 [HI [13 [H]. However, our algorithm is much faster 
than those methods in the literatures [13 [HI [IS [m for computing the asymmetric 
excited states. 


5. Concluding remarks. Different spatial discretizations including the finite 
difference method, sine pesudospectral and Eourier pseudospectral methods were 
adopted to discretize the energy functional and constraint for computing the ground 
state of Bose-Einstein condensation (BEC). Then the original infinitely dimensional 
constrained minimization problem was reduced to a finite dimensional minimization 
problem with a spherical constraint. A regularized Newton method was proposed by 
using a feasible gradient type method as an initial approximation and solving a stan¬ 
dard trust-region subproblem obtained from approximating the energy functional by 
its second-order Taylor expansion with a regularized term at each Newton iteration as 
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Fig. 6. Contour plots of the ground state ( f)g (a), asymmetrie exeited state in the x-direction 
010 (b), (c) asymmetrie exeited state in the y-direetion 0oi (c), an d as ymmetrie exeited state in 
both X- and y- directions 0ii (d) of a BEC in 2D with the potential jd = 500 in ^4-4- 



well as adopting a cascadic multigrid technique for selecting initial data. The conver¬ 
gence of the method was established by the standard optimization theory. Extensive 
numerical examples of non-rotating BEC in ID and 3D and rotating BEC in 2D 
with different trapping potentials and parameter regimes demonstrated the efficiency 
and accuracy as well as robustness of our method. Comparison to existing numerical 
methods in the literatures showed that our numerical method is significantly faster 
than those methods proposed in the literatures for computing ground states of BEC. 
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